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The heterogeneous force networks in static granular media — formed from contact forces between 
grains and spanning from boundary to boundary in the packing — are distinguished from other 
network structures in that they must satisfy constraints of mechanical equilibrium on every ver- 
tex/grain. Here we study the statistics of ensembles of hyperstatic frictionless force networks, which 
are composed of more forces than can be determined uniquely from force balance. Hyperstatic force 
networks possess degrees of freedom that rearrange one balanced network into another. We con- 
struct these rearrangements, count them, identify their elementary building blocks, and show that 
in two dimensions they are related via duality to so-called floppy modes, which play an important 
role in many other aspects of granular physics. We demonstrate that the number of rearrangements 
governs the macroscopic statistical properties of the ensemble, in particular the macroscopic flucu- 
tations of stress, which scale with distance to the isostatic point. We then show that a maximimum 
entropy postulate allows one to quantitatively capture many features of the microscopic statistics. 
Boundaries are shown to influence the statistics strongly: the probability distribution of large forces 
can have a qualitatively different form on the boundary and in the bulk. Finally, we consider the role 
of spatial correlations and dimension. All predictions are tested against highly accurate numerical 
simulations of the ensemble, performed using umbrella sampling. 
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Granular systems are athermal and dissipative: an un- 
driven system will eventually reach a static, mechanically 
stable state and remain there. By repeatedly applying 
the same preparation protocol, whether numerically or 
experimentally, it is possible to build up an ensemble of 
static granular packings in which each element of the en- 
semble is a final state of the preparation protocol. Such 
an ensemble is very different from an equilibrium ensem- 
ble: not only is there no thermal equilibrium, there are 
no dynamics at all. 

Micromechanical models like soft spheres interacting 
via repulsive contact forces are an attractive way to study 
granular materials numerically; see, e.g., Ref. PQ and ref- 
erences therein. Even at this level of abstraction, how- 
ever, a theoretical description of the statistical or me- 
chanical properties of disordered packings is daunting. 
In recent years, a model system called the force network 
ensemble (FNE) [5] has proven to be an extremely use- 
ful tool for studying the disordered stress states of static 
granular packings. In its simplicity the FNE affords the- 
oretical traction and permits highly accurate simulation, 
and yet it is detailed enough to reproduce many of the 
statistical and mechanical properties of numerical and 
experimental packings [3 9 ; Ref. [10] provides a review. 

The force network ensemble is built on the observation 
that packings of noncohesive frictionless disks or spheres 
at finite pressure are hyperstatic: they possess a number 
of contacts in excess of that which would uniquely de- 
termine the contact forces from mechanical equilibrium. 
This means that for a single packing geometry there exist 
many different configurations of forces that satisfy force 
balance on each grain. This indeterminacy can, in prin- 
ciple, be lifted by specifying a contact force law, from 



which the forces can be determined given the grain posi- 
tions. The conceit of the FNE, however, is not to specify 
this information, and instead to exploit the force inde- 
terminacy to practical advantage. Because deformations 
are so small in packings of hard but not perfectly rigid 
grains, there is a strong separation between the grain 
scale and the contact scale. The idea is that, by aver- 
aging over all balanced force networks on a single frozen 
contact network, one captures fluctuations in the stresses 
that would also result from rearrangements of the grain 
positions. In the FNE, then, grains do not rearrange but 
forces do (Fig. [I]) . 

Here we describe the stress statistics of the force net- 
work ensemble within a statistical mechanics framework. 
This work represents an elaboration and significant ex- 
pansion on recent results, which demonstrated that the 
statistics of local stresses within the FNE can be de- 
scribed using a maximum entropy principle [5]. The no- 
tion of an ensemble description of static granular matter 
dates to Edwards, who proposed an ensemble of packings 
characterized by like boundary conditions |llj . Though 
conceptually appealing, the Edwards ensemble is difficult 
to probe theoretically. In this spirit the FNE can be seen 
as a restricted but more accessible version of Edwards' 
ensemble. In ensuing years, a number of authors have 
proposed granular ensembles in which stress plays a role 
similar to that of energy in an equilibrium ensemble |12l - 
119] . and we shall see that the FNE naturally lends itself 
to such an approach. 

This paper is divided in two main sections. Section [T] 
develops the FNE in greater detail, with particular em- 
phasis on the force rearrangements that transform one 
force network into another. Force rearrangements can 
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be constructed "by hand" by considering the local force 
balance constraints on grains in a static packing. Their 
onstruction allows us to motivate the ensemble differently 
than previously [51 [TO] , in a manner that anticipates the 
theoretical framework to follow. Whereas in prior work a 
constraint on the global stress was explicitly imposed as 
a separate constraint, we demonstrate that the ensemble, 
defined in terms of force rearrangements, can be built up 
in a manner that obviates this step. It turns out that 
the force rearrangements act in such a way as to leave 
the global stress unchanged, without imposing an exter- 
nal constraint. This leads naturally to the idea of the 
FNE as an analog of the microcanonical ensemble: it is 
a system with "dynamics" , here realized by the force re- 
arrangements, that leave a quantity invariant. We will 
demonstrate that in the FNE that quantity is stress or a 
related invariant, rather than energy. 

Given its similarity to the equilibrium microcanonical 
ensemble, it is natural to anticipate that the ensemble 
dynamics maximize an entropy. In Section [TT] we write 
down this entropy and develop a framework that allows 
for description of the statistics of macroscopic stresses 
in the ensemble: We derive the FNE equation of state 
and an expression for macroscopic pressure fluctuations 
in the canonical FNE. We then turn to the statistics of 
local measures of the stress, particularly pressure p at 
the grain scale; our focus is on the form of the pressure 
probability distribution function of p for asymptotically 
small and large p. We also treat spatial correlations and 



dimensions d > 3, respectively. Finally, Section III gives 
a discussion and future outlook. 



I. FORCE REARRANGEMENTS AND THE FNE 

Energy ceases to be an invariant of the dynamics in a 
dissipative system. In ensembles of static granular pack- 
ings, moreover, we explicitly decline to consider the dy- 
namical states that the system traverses when it passes 
from one static state to another. In this context the force 
network ensemble can be understood as an intermediate 
case. As described above, each force network in the en- 
semble is meant to represent the stress state of a static 
packing. At the same time, we shall see that the ensem- 
ble does have a discrete dynamics, at least in the sense 
of Monte Carlo rearrangements of the forces. In this 
section we will construct these force rearrangements for 
arbitrary disk packings and then use them to give a more 
precise definition of the force network ensemble. We then 
demonstrate that the Monte Carlo dynamics of the FNE 
leave two related quantities invariant. These quantities 
are then natural candidates on which to base a statistical 
description; this is the topic of Section [TT] 

To begin, we illustrate the concept of force rearrange- 
ments in a particularly simple contact network, the tri- 
angular lattice (Fig. [TJ. Although the contact network 
is ordered, it admits many different force balanced con- 
figurations of the contact forces between the grains, the 
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FIG. 1: A single hyperstatic contact network, such as the pe- 
riodic frictionless triangular lattice, admits many force net- 
works in which forces balance on every grain; here three are 
shown. Lines represent contact forces, with their thickness 
proportional to the force magnitude. Though the contact net- 
work is ordered, typical force networks are disordered. 



overwhelming majority of which are disordered. This is 
because the triangular lattice is strongly hyperstatic. 

Hyperstaticity can be quantified in terms of the mean 
coordination number of a packing z; in the triangular 
lattice every grain has six contacts, and z = 6. A generic 
disk or sphere packing must have a minimum coordina- 
tion z; so to satisfy force balance on each grain. To see 
this, one can think of the individual contact forces as 
degrees of freedom. There must be enough degrees of 
freedom to satisfy all the constraints of mechanical equi- 
librium. Let us consider frictionless spheres, for which all 
contact forces have just one component, a normal force, 
and hence torque balance is satisfied automatically. In 
a frictionless packing of N grains in d dimensions, these 
constraints are the equations of force balance, of which 
there are Nr, = dN. There is one contact force for each 
of the N c — \zN contacts. A system is called isostatic 
when there are just enough forces to satisfy the con- 
straints. Requiring Nr, = N c yields z- lso ~ 2d. The cal- 
culation can be repeated for frictional spheres, in which 
case there are both additional constraints due to torque 
balance and additional degrees of freedom due to the tan- 
gential force components at each contact. The bulk of 
this work is dedicated to frictionless packings, although 
many results have straightforward generalizations to the 
frictional case. 



A. The triangular lattice 

The frictionless triangular lattice has an excess coordi- 
nation Az :~ z — z- lso = 2. Because each contact is shared 
by two grains, this means that there is \Az — 1 degree 
of freedom per grain. Because the triangular lattice is or- 
dered, it is possible to construct this degree of freedom by 
inspection. It is shown in Fig. [2j and known as a "wheel 
move" due to its appearance [4j. A wheel move is a re- 
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wheel moves to fo does not change the ensemble defined 
according to Eq. ([I]) [51]. We show in Section IF that 



FIG. 2: Monte Carlo or "wheel" move in the frictionless tri- 
angular lattice [4]. Force is added to spoke contacts (solid 
bars) and subtracted from six rim contacts (open bars), or 
vice versa. The net vector force on each of the seven partici- 
pating grains is unchanged. 



arrangement of 12 forces: six "spoke forces", which are 
changes to the contact forces on a central grain, and six 
"rim forces" , which change the contact forces between 
the nearest neighbor grains. The spoke and rim forces 
have equal changes in magnitude but opposite sign. By 
design, all seven grains are in force balance. Because the 
equations of force balance are linear, a wheel move can 
be added to any balanced force network on the triangular 
lattice, such as the ones in Fig. [T| without violating force 
balance. 

By labeling each contact force /y = —fji by the grains 
it acts on, we can construct a vector f = {fij} containing 
all the unique forces in the force network. One possible 
force balanced configuration in the frictionless triangular 
lattice is the force network fo in which every force has 
the same magnitude /. In the same notation, we label 
the N wheel moves, one for each grain, as {Sfk}, where 
the index indicates the grain on which they are centered. 

We now define the the (isotropic) force network ensem- 
ble on the frictionless triangular lattice to be all noncohc- 
sive force networks {fk} that can be reached by applying 
any combination of wheel moves to the force network fo. 
That is, a force network in the FNE can be expressed 



f = fo 



JV-l 
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w k Sf k . 



(1) 



The weights {uik} then serve as coordinates of the force 
network f in a high-dimensional space. Only N — 1 wheel 
moves are required for a linearly independent set [4], 
hence the upper bound in the sum. To see this, consider 
the application of every wheel move with equal weight 
- the changes to each force sum to zero. Therefore the 
{wk} are only defined modulo their sum or, equivalently, 
one of the rearrangements is dictated by the sum of the 
others. 

The restriction to noncohesive forces means that the 
weights {wk} must be chosen so that the normal force 
component (f n )ij > for all contacts; recall that in 
frictionless sphere packings all forces are normal forces. 
From the linearity of the constraints it follows that the 
space of force networks is a convex polytope [HI [20l [21] . 
Therefore the choice of fo is not essential: replacing f 
by any force network f^ that can be reached by applying 



the reference network f serves to select the stress tensor 
a. 

Dynamics in the FNE can be visualized as a random 
walk in the space of force networks. As their name antic- 
ipates, the wheel moves serve as Monte Carlo moves. To 
sample the space of force networks, Monte Carlo moves 
are randomly selected and used to update the current 
force network . The size of the move is uniformly selected 
from the interval of possible step sizes, determined by the 
positivity constraint on the affected contact forces |52) . 
In this way, for sufficiently long runtimes, the space of 
force networks is sampled with flat measure; see Refs. [8] 
and [4] for details. 

We now turn to constructing force rearrangements in 
disordered contact networks. With the exception of Sec- 
tion HE this work is concerned primarily with two- 
dimensional systems. We therefore illustrate the con- 
struction of force rearrangements in two-dimensional con- 
tact networks, where we have access to a geometric con- 
struct called a Maxwell-Cremona diagram or reciprocal 
tiling [22] . The reciprocal tiling is a helpful tool and, 
furthermore, allows us to make connections to the con- 
cept of floppy modes, which play an important role in 
other aspects of granular physics [T[ [23J HI] • We there- 
fore describe the reciprocal tiling first before addressing 
disordered force rearrangements. We emphasize, how- 
ever, that force rearrangements can also be constructed 
in higher dimensions [35] . 



B. Maxwell-Cremona diagrams 

Contact forces in a packing may be used to construct 
a Maxwell- Cremona diagram, in which pairs of action- 
reaction forces between grains are mapped to a tiling of 
the plane. An edges in the tiling has a length propor- 
tional to the magnitude of the corresponding force, and 
its orientation is perpendicular to the vector force. This 
is most easily seen graphically (Fig.[3^,-c): the boundary 
of a tile is constructed by rotating the vector forces acting 
on a grain by tt/2 and placing them end to end in a right- 
hand fashion. Because the boundary is the vector sum 
of the contact forces acting on the grain, the boundary is 
closed (a polygon) whenever the grain is in force balance 
and not subject to body forces. Though the tiling can be 
generalized to incorporate body forces |26| , they will not 
be considered here. 

By Newton's third law, tiles of contacting grains have 
faces of like length and orientation. Hence tiles may be 
placed next to each other seamlessly, and the Maxwell- 
Cremona diagram is built up tile-by-tile in this fash- 
ion. For a packing in static force balance subject 
to imposed forces at the boundary, the corresponding 
Maxwell-Cremona diagram has a closed boundary and 
no internal gaps, as illustrated in Fig. [3}i and e. For 
a periodic packing, the tiling is also periodic. The re- 
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FIG. 3: Constructing a tile in the Maxwell-Cremona diagram, 
or reciprocal tiling, (a) A disk with four contacts. Vector con- 
tact forces imposed by the neighboring disks are indicated by 
arrows. Note that the forces need not be frictionless, i.e. need 
not be parallel to the segment connecting the contacting disks' 
centers, (b) Each contact force is rotated by ^ and drawn 
end-to-end, proceeding around the grain in a right-hand fash- 
ion, (c) The polygon enclosed by the vectors is the grain's 
corresponding tile in the tiling. The corners of the tile are 
vertices in the Maxwell-Cremona diagram. (d,e) Due to New- 
ton's third law, tiles corresponding to contacting grains can 
be placed flush against one another. The result is a tiling or 
tesselation. To determine the area of the tiling, it suffices to 
know the boundary forces on the packing; alternatively, one 
can sum the areas of individual tiles. 

ciprocal space coordinates of the diagram's vertices are, 
after rotation by — 7r/2, the void forces of Satake [27] and 
equivalently the loop forces of Ball and Blumenfeld [5H] . 

The reciprocal tiling exists as a consequence of static 
force balance in the packing. The construction makes no 
assumptions regarding the presence or absence of tan- 
gential or tensile forces, and torque balance is not a nec- 
essary condition for its existence. Nevertheless, all the 
force networks we study here are noncohesive and (triv- 
ially) torque balanced. 

For later convenience let us now calculate the mean 
coordination number in a periodic reciprocal tiling. For 
each grain, contact, and void in the packing there is a 
corresponding tile, edge, and vertex, respectively, in the 
tiling, and the topology of the tiling is the same as that of 
the dual graph of the contact network [53]. Euler's rela- 
tion for a periodic network relates the number of grains, 
contacts, and voids in a packing: 

N - N c + N = , (2) 

where N is the number of voids, or equivalently the num- 
ber of vertices in the reciprocal tiling. The mean coor- 
dination number in the tiling is z = 2N C /N. Combined 
with Eq. ([2]) this gives 

7y = - + - , 3 
2 z z 



FIG. 4: Two vertices in a reciprocal tiling, (a) Floppy mo- 
tion displacements (arrows) may rotate the edge joining the 
vertices (thick dashed line) but not change its length, (b) Ro- 
tating the displacements from (a) by n/2 produces a motion 
that changes the edge's length, but not its orientation. 

where we have written the relation in a form that em- 
phasizes the symmetry between z and z. In particular, 
whenever z > z- lso , z < z; so . 

C. Floppy motions and force rearrangements 

We now show that there is an intimate connection be- 
tween the force rearrangements in frictionless packings 
and floppy modes — non-rigid body motions of a collec- 
tion of particles that do not change the inter-particle dis- 
tances. The connection enters via the Maxwell-Cremona 
diagram. A motion of the tiling vertices that does not 
change the distance between vertices is a floppy mode. 
Let us label the coordinates of vertices in the reciprocal 
tiling as {hi}, i = 1 . . . N. If there is a floppy mode in the 
tiling, it must be that for each vertex with coordinate hi 
connected by an edge to a vertex at hj , their motions Shi 
and Shj are such that 

(hi - hj) ■ (Shi - Shj) = 0. (4) 

That is, to leading order in the motions, the relative mo- 
tion of the two vertices can have no component along the 
original orientation of the edge between them; this is il- 
lustrated in Fig. [4^. The same consideration applies to 
every edge in the tiling. 

The key observation is as follows, and is illustrated in 
Fig. [4] A given floppy mode describes a set of vertex mo- 
tions {Shi}. We take these vertex motions, rotate each 
of them by tv/2, and label the new motions {8h^}. By 
Eq. Q, the new relative motions are such that, to lead- 
ing order, the edge between any pair of connected ver- 
tices does not rotate; instead it translates and changes its 
length. This is most easily seen in Fig. [4] Recall, how- 
ever, that the edges in the reciprocal tiling are the forces 
in the original packing, and that a force rearrangement 
in a frictionless packing must change the magnitude but 
not the orientation of a contact force. Moreover, because 
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FIG. 5: A portion of a soft disk packing; lines indicate edges 
in the Delaunay triangulation of the disk centers. Edges with 
thick lines are also a part of the packing's contact network. 




FIG. 6: Disordered wheel moves, the building blocks of fric- 
tionless force rearrangements, (a) In a triangulation, every 
grain/ vertex is surrounded by a closed loop formed by its 
nearest neighbors. We define angles 8 and ij) for each triangle 
around the grain, (b) The grain's reciprocal tile has just one 
dangling edge from each vertex. A local floppy motion can be 
implemented by moving each vertex of the tile perpendicular 
to the dangling edge with a magnitude Shi determined in the 
text, (c) The disordered wheel move can be constructed by 
rotating the floppy motions {Shi} by ^. Moves in contact 
networks that are not triangulations can be constructed from 
linear combinations of these wheel moves. 

the motions are performed in the reciprocal tiling, which 
only exists if the packing is in force balance, they must 
respect force balance. Therefore floppy modes in the re- 
ciprocal tiling can be mapped to force rearrangements in 
the packing simply by rotating their motions by n/2. We 
will use this fact to construct a set of "disordered wheel 
moves" . 



D. Rearrangements in triangulations 

We will construct force rearrangements in disordered 
contact networks in two steps. In the first step we con- 
sider contact networks that are periodic triangulations of 
the plane. There it is possible to identify straightforward 
generalizations of the wheel moves in the triangular lat- 
tice. In the second step we construct linear combinations 
of these local rearrangements to form force rearrange- 
ments on contact networks that are not triangulations. 

Given the grain positions in a disk packing, the pack- 
ing's Delaunay triangulation can be constructed [30] ; see 
Fig. [5] for an example. A Delaunay triangulation con- 
nects disks that are geometric neighbors, regardless of 
whether they are actually in contact. In systems with 
large polydispersity there may be physical contacts that 
are not Delaunay edges; it is possible to correct for this, 



but here we assume that all physical contacts are indeed 
edges in the triangulation. An important feature of tri- 
angulations for present purposes is the following: each 
vertex is enclosed by a loop formed from edges connect- 
ing the neighbors of the vertex, i.e. the vertices to which 
the central vertex is connected by an edge. This is de- 
picted in Fig. [6^,. The idea is to first treat all edges in the 
triangulation as if they can carry force. We then iden- 
tify a set of floppy modes centered on each tile in the 
tiling; the floppy modes then give the disordered analog 
of a wheel move by the above prescription. In general, 
the wheel moves will make changes to the force on edges 
that are part of the Delaunay triangulation but not the 
contact network. To correct for this, we later construct 
linear combinations of the wheel moves that do not alter 
the force on these "deleted edges" . 

Fig. [6Jd depicts the reciprocal tile corresponding to the 
grain in Fig. [6^,. Because the neighboring grains form a 
closed loop about the central grain, the central grain's 
corresponding tile has just one edge extending from each 
of its vertices. These "dangling edges" allow a floppy 
mode to be constructed by displacing the vertices of just 
a single tile. For a tile with Zj vertices, each with coor- 
dinate hj in the reciprocal space, we prescribe displace- 
ments {5hj}, j — 1 . . . Zj. Each Sh must be orthogonal 
to the dangling edge at its vertex, in accordance with 
Eq. The magnitudes {Shj} of the vertex motions 

can be determined from the constraints that the lengths 
of the tile's edges remain unchanged to leading order in 
the motions. Although there are magnitudes and as 
many constraints, this is a one parameter family of mo- 
tions: the constraint equations are homogeneous in the 
magnitudes {Shi}, hence there is a degeneracy among 
them. 

To see the degeneracy, we apply Eq. Q to each edge 
around the boundary of the tile. This gives Shj sinOj = 
Sh,j + i sini/jj, where the angles 9 and ip are defined in 
Fig. [6] and the index j increments around the tile vertices 
in a righthand fashion, modulo Zj. Taking the product 
of this relation applied to each edge of the tile, one finds 

z z 

Y[ sin 9i = Y\ sin ^ . (5) 

i=l »=1 

Note that Eq. ^ involves only angles; the magnitudes 
of the vertex displacements have dropped out. In fact, 
the same relation arises by repeatedly applying the law of 
sines to the triangles surrounding the central node in the 
original triangulation (Fig. [6^,). This is a consequence 
of the fact that edges between nearest neighbors of the 
central grain form a closed loop. It means that Eq. ^ 
is satisfied automatically by geometry, and only z, — 1 
of the Zj apparent constraints on the magnitudes {Shi} 
are independent. The resulting one parameter family of 
motions is a floppy mode in the Maxwell-Cremona dia- 
gram; it is, by construction, localized to the vertices of a 
single tile. By following the prescription {Shi} {Sh^} 
we have thus constructed a localized force rearrangement 
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(Fig.^). There is one per grain; we call them disordered 
wheel moves and label them {8fJ° ca1 }, k — 1 . . . N. 

The force network ensemble on a triangulation can thus 
be defined to comprise all force networks of the form 



Each such rearrangement can be written 



f = fo 



JV-l 

£ 



local 



(6) 



The weights of the disordered wheel moves {wk} again 
serve as coordinates of the state in the space of force 
networks. The sum runs to N — 1 because, as in the 
triangular lattice, there is a sum rule on the disordered 
wheel move weights (see Appendix A). The particular 
solution fo can be identified via simulated annealing [31 
or linear programming [4]. 

The reciprocal tiling is similarly useful in constructing 
the counterpart of disordered wheel moves in frictional 
systems. This procedure is described in Appendix B. 



E. Rearrangements in general disordered contact 
networks 

Let us now consider rearrangements in periodic contact 
networks that are not triangulations. Euler's relation, 
Eq. ([2]), guarantees that a packing's Delaunay triangula- 
tion has mean coordination number z tr ; = 6. By com- 
parison, the physical contact network with mean coordi- 
nation number z has only N c = \zN contacts, so that 
there are = 3N — N c deleted edges in the triangula- 
tion, i.e. edges in the triangulation that cannot transmit 
force because they do not correspond to contacts in the 
packing. Therefore we must impose constraints on 
the disordered wheel moves from the triangulation: they 
must always act in combinations chosen so as to neither 
add nor subtract force on deleted edges. We now give a 
prescription for constructing these combinations. 

Define F to be the N - 1 x \zN matrix with 5f l k ocaX , 
k = 1 ... N — 1, as its fc th column. Likewise define D to 
be a \zN x 7V<j matrix such that each column has a unit 
entry for a unique deleted edge, and all other elements 
zero. Then, for w = {wk}, the matrix operation (DF) w 
returns the effect of a superposition of local moves with 
weights w on the deleted edges. As the deleted edges 
cannot carry force, this must be null: 



(DF)w = 0. 



(7) 



A set of vectors {w„}, n = 1 . . . N w , spanning the null 
space of DF gives a basis set of linear superpositions of 
the disordered wheel moves {(5f^ oca1 } having null effect on 
the deleted edges. These are our force rearrangements. 
Assuming none of the deleted edges are redundant, which 
is generally true for a disordered network, each deleted 
edge reduces the number of force rearrangements by one 
from the N— 1 independent disordered wheel moves avail- 
able in a triangulation: N w = N — JVj — 1, or 
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N-l 

k=l 



(9) 



where [w„]fc is the k th element of w„. The label empha- 
sizes that these rearrangements are spatially extended, 
being composed of multiple disordered wheel moves. 

In analogy to Eq. ([I]), we can now define the FNE on 
a disordered contact network as the set of force networks 
expressible as 



f = fo 



iV w 

•E 

71=1 



r 8f cx 

c 7i tJL n 



(10) 



-AzN- 1. 



(8) 



where the iV w coefficients {c n } now serve as coordinates 
of the force network. Typically the protocol used to gen- 
erate the contact network - here we employ molecular 
dynamics simulations - also produces a force network 
that can be used as the particular solution fo. As in a 
triangulation, alternative choices for fo can be identified 
via simulated annealing or linear programming. 



1. Relation to the isostatic length 

Because force rearrangements are closely related to 
floppy modes, it is natural to expect a connection to the 
"isostatic length" I*. The isostatic length governs sev- 
eral mechanical properties of packings of soft frictionless 
particles [U [24J I3TH34] and is related to floppy modes via 
the "bond cutting" argument of Wyart and co-workers 
[251 |2H d2] ) which we briefly summarize. Consider an in- 
finite or periodic packing with mean coordination num- 
ber z = z iso + Az. z > z iso so that the packing has 
no floppy modes. Now imagine removing a cluster of 
size I from the packing, as in Fig. [3H; in so doing one 
cuts 0(£ d_1 ) "bonds", viz. contacts, on the boundary. 
The rigidity of the cluster is determined by a competi- 
tion between cut bonds, which remove constraints and 
therefore inhibit rigidity, and "excess bonds" (in excess 
of an isostatic packing) in the interior, which lend redun- 
dant rigidity to the cluster. If there are more cut bonds 
than excess bonds, of which there are 0(Az£ d ), the clus- 
ter will possess floppy modes. The marginal case occurs 
when the two populations balance, and this selects the 
isostatic length I* ~ 1/ Az. 

Let us now consider the Maxwell-Cremona diagram of 
the same packing. Its mean coordination number is be- 
low isostaticity, cf. Eq. ([3]), and so the diagram contains 
floppy modes. We again select a cluster of grains from the 
packing, only now we fix the 0(i' d_1 ) forces on its bound- 
ary. In the reciprocal tiling this isolates one portion of 
the tiling from the rest, see Fig. [3^. Fixing forces adds 
constraints, and if there are enough of them to overcome 
the shortfall of bonds in the interior (compared to an iso- 
static packing), the isolated region will be rigid. As the 
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(b) 




Fi 



FIG. 7: (a) A force network from the triangular lattice and 
(b) its Maxwell-Cremona diagram or reciprocal tiling. 



shortfall of bonds is 0(\Az\ N) ~ 0{Az£ d ) and the num- 
ber of fixed forces is a boundary term, the isolated region 
will be rigid when £ is sufficiently small. The marginal 
case is again £* ~ 1/Az. Recalling that floppy modes 
in the reciprocal tiling construct force rearrangements, 
we see that the isostatic length can be interpreted as the 
typical linear size of a grain cluster that supports just 
one internal force rearrangement. Indeed, it was recently 
shown in Ref. |35j that point-to-set correlation functions 
in the FNE, which measure the overlap between force 
networks in a finite size cluster, scale with I*. 



F. Additive Extensive Invariants 

We have now shown that in hyperstatic frictionless disk 
packings it is possible to construct a set of contact force 
rearrangements that transform one force balanced force 
network into another. The force network ensemble can 
then be defined as all noncohesive force networks that 
can be reached, starting from some initial force network 
fo, via these force rearrangements. We now show that 
the only salient feature of the initial force network is its 
stress tensor a. In so doing, we will identify two additive, 
extensive quantities that are invariant under the Monte 
Carlo dynamics of the force network ensemble. 

The reciprocal tiling was helpful in constructing force 
rearrangements; we use it again here to describe the 
stress tensor of a statistically homogeneous force net- 
work. Consider the periodic force network of Fig. [7^, and 
its reciprocal tiling in Fig. [7]3. The packing has orthogo- 
nal primitive vectors L\ = L\e\ and L2 — L2&2, and any 
surface normal to e+ experiences a net compressive force 
Fi = aLi. Without loss of generality we choose e\ and 
e.2 to align with the principal stress directions, so that 



o-ii 
CT22 

<5l2 



Li 
F2 
L 2 
0. 



(11) 



At the same time, by periodicity the tiling must have 
primitive vectors F\ and F 2 . Because the contact net- 



work is fixed in the FNE, the "size and shape" of the 
tiling directly encodes the stress tensor via Eq. (111. In 



the following sections, the tiling's area A will play an 
important role. Note that 



A = F X F2 
= (deta)V, 



(12) 



and is therefore extensive. From the construction of the 
tiling, A is manifestly additive: 



N 



(13) 



A = a,, 

where a iy the area of the tile corresponding to grain i, is 
1 p 

ai = ^e 3 -22gjXg j+1 . (14) 
3=1 

The sum runs over the Zi contacts of grain i ordered in a 
righthand fashion around the grain, and indices are taken 
modulo Zi. The vector cjj+i ■= X}i=i fjk- Note g\ = 0. 
The unit vector £3 points out of the plane in a sense such 
that a, is positive when all forces are compressive. 
Eq. (Ill has important implications for the force net- 



work ensemble. To change the stress tensor or tiling area, 
a force rearrangement must change the primitive vectors 
of the reciprocal tiling's unit cell. A wheel move, disor- 
dered or not, cannot change the primitive vectors F\ and 
F 2 . To see this, consider the action of a wheel move in 
the tiling, as in Fig. [6];: the move simply shuffles area 
among a small number of tiles, leaving the unit cell of 
the tiling unchanged. Therefore the stress tensor a, or 
equivalently the extensive stress S :— aV, is a topolog- 
ical invariant of the FNE, as is the tiling area A. This 
observation applies equally to general disordered contact 
networks, because there force rearrangements are linear 
superpositions of disordered wheel moves. 

The force network ensemble therefore bears strong sim- 
ilarity to the microcanonical ensemble. Just as energy is 
an additive, extensive invariant of the dynamics in an 
equilibrium system, the extensive stress S is an additive, 
extensive invariant of the Monte Carlo dynamics of the 
FNE. Moreover, because "dynamics" in the FNE are per- 
formed by a random walk in the space of force networks, 
force networks are sampled with equal a priori probabil- 
ity. This is again reminiscent of the equilibrium micro- 
canonical ensemble. In the following, both for simplicity 
and to reinforce the analogy, we restrict our attention to 
ensembles of isotropic force networks, so that the scalar 
"extensive pressure" V — ^Tr6> fully specifies the stress, 



S = V1. 

The extensive pressure is also additive, 



N 



(15) 



(16) 
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where 



1 N 

y 

3 = 1 



fij 



(17) 



and fij is the force on grain i applied by grain j (nonzero 
only when i and j are in contact) and fij is the vector 
from the center of j to i. We pursue the analogy to 
equilibrium ensembles further in the following Section. 

Note that in previous work, the FNE has been defined 
as the flatly sampled ensemble of noncohesive force net- 
works subject to local force balance constraints and an 
imposed stress tensor a [2J [3] . The definition of the force 
network ensemble we have offered is ultimately equiva- 
lent (see Appendix A), but instead of directly imposing 
a conserved stress, we arrive at it naturally through con- 
sideration of local rearrangements consistent with force 
balance. From this perspective, the comparison to the 
microcanonical ensemble, and the statistical mechanics 
analogy developed in the following section, is more apt. 

If the usual definition of the FNE is a microcanonical 
one, can one pass to a canonical FNE? This question was 
considered in detail in Ref. [36]; the answer is yes. In a 
canonical ensemble the extensive pressure V is not invari- 
ant but is allowed to fluctuate. According to the above 
discussion, this cannot be achieved with superpositions of 
local force rearrangements. A canonical FNE, therefore, 
must admit additional force rearrangements that change 
the global stress of a force network. Still restricting our- 
selves to isotropic stress states, just one additional rear- 
rangement is needed; one possibility is the rescaling of all 
the forces in a force network by a single scalar parameter. 
In the reciprocal tiling this corresponds to an affine in- 
flation or contraction of the entire tiling. In a canonical 
ensemble, of course, not all force networks receive equal 
weight; we derive the equivalent of the Boltzmann factor 
in the following Section. 



II. STATISTICS IN THE FNE 

In the previous section we demonstrated that the force 
network ensemble can be built up from force rearrange- 
ments that respect the constraints of local force balance. 
We now turn to a study of stress statistics in the ensem- 
ble we have constructed. After writing down and maxi- 
mizing entropy in the FNE, we consider pressure in the 
canonical ensemble, deriving the equation of state and 
the scaling of pressure fluctuations. We then consider 
the statistics of stress at the grain scale in the form of 
the local pressure probability distribution P(p). 

Because of its similarities to equilibrium ensembles, it 
is natural to describe the force network ensemble within 
a statistical mechanics framework. By definition - or al- 
ternatively, the Monte Carlo dynamics described above 
guarantee that - every force network in the FNE is sam- 
pled with equal a priori probability, i.e. with a flat mea- 
sure. Here we will write down an entropy, postulate that 



it is maximized, and show that it correctly reproduces 
equal a priori sampling in the microcanonical FNE. We 
then follow the same approach to generate a canonical 
FNE. This approach is in direct analogy to the standard 
textbook treatment of equilibrium statistical mechanics 
[37], though we spell out the steps for completeness. 

Because both V and A are additive extensive invari- 
ants of the dynamics, and because they bear a one-to-one 
correspondence in noncohesive isotropic force networks, 



see Eq. (12), one can employ either as the analog of en- 



ergy in an equilibrium ensemble. This point is discussed 
at length in Ref. [3B]. Here we will use the extensive 
pressure V for greater similarity to other approaches in 
the literature [T2HI9] . 

If a force network f occurs with (normalized) frequency 
B({), there is an associated entropy 



S[B] 



df G(f) [5(f) In B({)} 



(18) 



The function G(f ) restricts the integral to "valid" states, 
and is defined to be unity when (i) f is force balanced, (ii) 
its contact forces are noncohesive and (Hi) the force net- 
work is isotropic, i.e. an = CT22 and CT12 = 0. G(f) = 
otherwise. The integral may be further restricted de- 
pending on the ensemble under consideration. We pos- 
tulate that S[B] is maximized subject to certain con- 
straints. 

Let us first consider the microcanonical ensemble, in 
which the relevant constraint is that of normalization: 



1 



df G{{)B(f). 



(19) 



V(f)=V 



The integral is restricted to force networks with extensive 
pressure Vq. The entropy S[B] can be maximized subject 
to Eq. (|19|) by introducing the Lagrange multiplier \ and 
writing 



= 6 1 dfG(f) 

>V(f)=V 



ln£(f)+x(l) 



B(f). (20) 



Eqs. (|19j) and (20) may be solved for x an d B(i); one 
finds B[f) = 1/0^0)1 where 



n(V ) = J diG(f)6(P(t)-P Q ) 



(21) 



Reassuringly, we find that valid force networks in the 
microcanonical FNE receive equal statistical weight. 

An extensive discussion of the microcanical FNE can 
be found in Ref. [3J. In the present work, we will find 
it convenient to adopt the perspective of a canonical en- 
semble, in which the extensive pressures V(i) is allowed 
to fluctuate; i.e. Eq. ( 19 ) is replaced by 



1 



dfG(f)5(f). (22) 
At the same time, a constraint on the average extensive 



pressure (V) is imposed, 

(V) = J df G(f)-P(f)5(f). 



(23) 
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We demonstrate below that the microcanonical and 
canonical ensembles are equivalent in the usual way. 

Maximizing entropy with the Lagrange multipliers £ 
and a, now subject to the constraints of Eqs. (22) and 
(23), we have 



= 6 df G(f) 



•lnB(f) + C(l) + aP(f) B(f). (24) 

The extremal distribution is 

B{{) = Z- 1 exp(-a'P(f)), (25) 

where Z := e^ 1 . The La grange multiplier a, the inverse 
of which was termed angoricity by Edwards |18j , is chosen 
to satisfy 



(V) 



d_ 

Oa 



\nZ. 



(26) 



and the partition function Z enforces normalization of 



Z{a)= J dfG(f) exp(-aP(f)). 



(27) 



Note that, in perfect analogy to the equilibrium case, 
force networks in the canonical FNE are weighted by a 
"Boltzmann factor" exp(— aP). 



A. Macroscopic quantities 

It is now straightforward to consider the statistics of 
the extensive pressure V, including an equation of state 
relating the intensive parameter a to (V). 

Our starting point is the partition function Z of 
Eq. (27), which may be re-expressed in terms of £l(V), 
the density of states with extensive pressure V: 



Z 



dVn(V) exp(-aV) . 



(28) 



fi('P) has a straightforward geometric interpretation: it 
is the content of the space of valid force networks with 
extensive pressure V. Recall that the space is a con- 
vex polytope in N w dimensions. V sets the typical "di- 
ameter" or linear dimensi on o f the polyt ope, so that 
tt(P) oc V Nv . From Eqs. ((26]), (|8j), and ([2l{, it then 
follows that 



a(V) = ^Az N + 0(1) 



(29) 



Up to corrections that vanish in the thermodynamic 
limit, Eq. (29) is the equation of state of the FNE. It 
simply states that a -1 selects the natural pressure scale 
in the ensemble. This was forseeable: having discarded 
the force law, the FNE does not possess an intrinsic force 
scale, leaving no other scale with which a -1 could com- 
pete. Eq. (29) is confirmed numerically in Fig. [8^,. 




60 

o 



3 4 

log 10 AzN 




log 10 &zN 



FIG. 8: (a) Equation of state for the FNE computed in the 
canonical ensemble for a = 0.05, 0.1, and 0.2. Solid curves are 
Eq. Q29n . For each data point a contact network was prepared 
by molecular dynamic simulation of a packing of N particles, 
with N ranging from 250 to 2000; the resulting and mean co- 
ordination numbers range from z = 4.25 to 6.00. Whenever 
necessary rattlers have been removed, (b) The relative fluctu- 
ations A 2 = {{SV) 2 )/{V) 2 for the same data as in (b) collapse 
when plotted against Az N, independent of a, as predicted by 

Eq. Jin |. 



The extensive pressure fluctuations can be calculated 
similarly, 

POO 

({6V) 2 ) = Z- 1 dVCl(V)(V- (V)) 2 exp(-aV) 
Jo 

= N w /a 2 . (30) 
Hence the relative fluctuations of pressure are 



,,_((SV) 2 )_ 2 



AzN 



(31) 



Note that the pressure fluctuations are governed by iV w , 
the dimension of the space of valid force networks. The 
equivalence of canonical and microcanonical ensembles 
follows from the 1/N scaling of the relative fluctuations, 
which vanish in the thermodynamic limit. Note the de- 
pendence on Az in Eq. (31), which ensures diverging 
relative fluctuations in the isostatic limit Az ! [53] . 
Eq. (31) is confirmed numerically in Fig. ^p, which plots 



the pressure fluctuations in simulations of the canonical 
FNE. The figure shows data for a range of system sizes N 
and coordination numbers z, all of which fall on the curve 



described by Eq. (31). In addition, for each (N,z) pair, 
three different values of a are plotted; these are difficult 
to see because, as predicted by Eq. (31), the fluctuations 
are independent of a. 

The scaling of pressure fluctuations can also be written 
A 2 = l/p w V. Therefore the pressure fluctuations are 
governed by the ratio of the linear system size C := V 1//d 



to the length scale £ w := 1/pw ~ l/Az 1 ^, namely 



C 



(32) 
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£ w sets the scale on which one finds fluctuations on the or- 
der of the mean pressure in periodic packings. We stress 
that £ w is different from, though closely related to, the 
isostatic length £* ~ 1/Az. The relative pressure fluc- 
tuations are not sensitive to the size of the force rear- 
rangements, but rather to their number. Their number 
scales as V/(£ w ) d rather than V/(£*) d . This is possible 
because there is significant spatial overlap between the 
force rearrangements, i.e. a typical contact participates 
in multiple rearrangements. 

Equivalence of the microcanonical and canonical force 
network ensemble can be expected only for periodic sys- 
tems of linear size C 3> £ w . In non-periodic systems 
the balance of boundary and excess bulk contacts again 
becomes relevant in constructing a canonical ensemble 
|36j - the system must at a minimum be large enough 
to support force rearrangements. The isostatic length 
I* ^> £ w as the isostatic limit Az I is approached, 
hence we anticipate the stricter requirement C £* 
for canonical-microcanonical equivalence in non-periodic 
systems. Nevertheless, the scaling of Eq. (32) must still 



within the microcanonical and canonical FNE, respec- 
tively. Here p\ is the pressure on the grain with index 
1 calculated via Eq. (fl7|). Formally the distribution is 



particular to the choice of grain on which the pressure 
Pi is assigned, unless the contact network is a Bravais 
lattice, though in practice little difference is found [3]. 
These two distributions of Eq. ( 33 1 converge to the same 



form in the thermodynamic limit, as demonstrated nu- 
merically in Ref . . Therefore to study the local stress 
statistics one may choose to work in whichever ensemble 
is convenient. 



Although Eq. ( 33 ) can be solved exactly for very small 



systems [31 11] , one must resort to asymptotic or approxi- 
mate methods in larger systems. We first show that P(jp) 
has a power law form for asymptotically small p, with an 
exponent dictated by local topology and force balance. 
We then develop an expression for the full distribution 
P(p); though the treatment is approximate, it is suffi- 
cient to capture quantitatively the pressure statistics in 
the frictionless triangular lattice. 



hold for asymptotically large non-periodic systems. 



B. Microscopic quantities 



Statistics of small pressures 



We now turn to the statistics of local stresses. One 
microscopic measure of the stress is the pressure p on 
an individual grain. Although the distribution of con- 
tact forces P(f) is widely studied, we will focus largely 
on the pressure distribution P(p), which has at least two 
advantages over P(f). First, p is slightly coarse-grained 
with respect to /, making it a more realistic target for 
the approximate expressions we develop in later sections. 
Second, it will prove to be convenient that p, being de- 
fined on the grain scale, enters at the same scale as the 
constraints of local force balance. 

The local pressure distribution is given by 

P„(p) = MVo)}- 1 J dfG(f)Wf)- Po)%i(f)-p) 
P a (p) = [Z(a)}- 1 J diG({)e- aV ^S( Pl ({)-p), (33) 



We now demonstrate that the local pressure distribu- 
tion on a grain with Zi contacts scales as P(pi) ~ p"* in 
the limit pi J, 0, where Vi = — d — 1 for frictionless 
spherical grains. In other words, for small pressures the 
statistics are governed by local topology and constraints. 

The notation f , used above, is a compact way to indi- 
cate the set of contact forces {fij} — {fij&ij}, where 
indicates the unit vector aligned with the contact from 
i to j. Here we employ a more explicit notation. In an 
enumeration of unique contact forces, each contact (i, j) 
need appear only once because fij = —fji by Newton's 
third law. In the following, however, we allow the prod- 
uct over contacts rifij) to include both and (j, i), 
and explicitly write down the constraint due to Newton's 
third law. 



PM = z- 1 / 

Jo 



dfye-^M* 6 (fij + fji) 



rn>2 \(i=m,j) 



where 



J(Pi,{fji}) 



U'J ,' '"" Kfij • fji) 
(ij) 



Jipuifji}) 



(34) 



The (5-functions serve to impose Newton's third law and vector force balance. The extensive pressure V in the 
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Boltzmann factor has been expressed in terms of the in- 
dividual contact forces. Terms involving the grain in- 
dexed 1, on which the pressure p\ is imposed, have been 
explicitly separated in the function J(pi, {fji})- 

Note that the only dependence on p\ enters through 
J(pi,{fji}). Secondly, the term e~ a ' Pl s» 1 for small p\. 
Further, grain 1 "interacts" with other grains only via 
Newton's third law. Thus, assuming interactions can be 
neglected in J(pi, {fji}) for asymptotically small pi, the 
pressure distribution scales as 




xS 



(35) 



The righthand side of Eq. ( 35 ) may be conceptualized 



in the following way. Given a grain with Z\ contacts, 
each force balanced configuration of the grain is a point 
in a state space in which the Z\ contact force magni- 
tudes are coordinate axes. Though the number of axes 
suggests Z\ dimensions, force balanced configurations oc- 
cupy a (z\ — (f)-dimensional volume because of the d force 
balance constraints on the grain. Single grain force con- 
figurations with a pressure p\ reside on a (z\ — d — 1)- 
dimensional "slice" through this volume. The integral of 
Eq. p5"| ) is the content of this slice. Because the state 
with all contact forces being zero (hence zero pressure 
Pi) is force balanced, and because all contact forces must 
be non-negative, the origin is a corner of the volume of 
balanced force configurations. Moving out of this corner, 
i.e. increasing px, the content of slices grows as p\ l d l ■ 
We thus have 



(36) 



for small pressures in frictionless packings. Employing 
the canonical ensemble did not play an important role; 
repeating the calculation in the microcanonical ensemble 
gives the same result. This is to be expected; in the ther- 
modynamic limit, the statistics of microscopic quantities 
should not depend on the choice of ensemble. 

We stress the implication of Eqs. (35) and (36 1: for 



weak interactions, the scaling of the pressure probabil- 
ity distribution can be inferred from degree of freedom 
counting — local contacts versus local force balance con- 
straints — and the flat measure on the ensemble. 

To test these results or equivalently the reasonable- 
ness of neglecting interactions in the above calculations, 
we perform numerical simulations of the FNE in a disor- 
dered packing. The conditional probability distribution 
P{p\z) is then sampled, i.e. the probability of obtaining 
a pressure p given that a grain has z contacts. Because 
the resulting curves are smoother, we plot the cumula- 
tive distribution C z (p) := J P dp'P(p'\z). Fig. [9] plots of 
log C z (p) versus (z — d) \ogp for a frictionless system; the 
linearity of the curves for small p confirms Eq. (j36h . 




(z-2) \og m p 

FIG. 9: Numerically sampled cumulative distributions 
C z (p) = Jq dp' P(p'\z). Statistics are sampled in the micro- 
canincal FNE on a periodic disordered packing in d — 2 di- 
mensions composed of 17 grains with 3 contacts, 233 grains 
with 4 contacts, 489 grains with 5 contacts, 279 grains with 
6 contacts, and 6 grains with 7 contacts. 



2. Statistics of large pressures 



The success of the approach in Section |II B 1| indicates 
that in the limit of asymptotically small pressures it is 
possible to adopt a "single grain picture" , i.e. to simplify 
calculations by neglecting interactions with neighboring 
grains, and still successfully predict local pressure statis- 
tics. Ref. [3] demonstrated that an approximate treat- 
ment of the statistics of local pressure in a single grain 
picture can be surprisingly accurate, even for p > (p). 
We now describe this approach. 

We demonstrated above that a maximum entropy pos- 
tulate correctly reproduces the appropriate equal a priori 
weighting of valid force networks in the microcanonical 
FNE, and used the same method to construct the canoni- 
cal FNE. All of the results derived in this manner have di- 
rect analogs in equilibrium statistical mechanics. We now 
employ the principle of maximum entropy more broadly. 
When all that is known about a system is that it must 
satisfy certain constraints, the "best guess" is that the 
system's state is described by a probability distribution 
that maximizes (information) entropy |38j . This state- 
ment reduces to the approach of Section |II A| if one im- 
poses all the relevant constraints on the system, includ- 
ing the constraints of local force balance on every grain. 
Though these local constraints are conceptually unprob- 
lematic and straightforward to implement in simulation, 
they render expressions like Eq. (34) difficult or impos- 
sible to evaluate. Because of these technical difficulties, 
we will seek to make approximations. In so doing, we 
take a simple lesson from information theory: the more 
information one incorporates (in the form of constraints 
on the system), the more accurate one can expect the 
predicted pressure distribution to be. 

We now perform a calculation in which entropy is max- 
imized subject to constraints on both (V) and (A). Wc 
have seen that, in the presence of local force balance, 
these two constraints are redundant. However, we will 
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also assume a single grain picture in which interactions 
with neighboring grains are not explicitly incorporated. 
In so doing, the mechanism by which the constraints on 
(V) and (A) are redundant is broken — there is no tiling 
unless every grain is in local force balance. In this con- 
text, imposing the constraint on (A), in addition to (P), 
reintroduces some of the information that was lost by ne- 
glecting interactions. In effect, it tells the central grain 
something about the consequences of force balance on 
all the other grains in the system. The surprise is that 
with this one additional piece of information it is possi- 
ble to describe local pressure distributions quantitatively. 
There is a price to be paid for this approach: because (V) 
and (^4) are not truly independent constraints, the La- 
grange multipliers a or 7 that we introduce to impose 
them cannot be associated with "true" intensive thermo- 
dynamic parameters in the thermodynamic limit. There- 
fore a and 7 within this approximate method should not 
be invested with any physical significance. 

We will treat the case of the frictionless triangular lat- 
tice. The main simplification comes from the ordered 
contact network; we saw above that the pressure distribu- 
tion P(p) depends on the local coordination number. In a 
Bravais lattice, of course, each grain has the same number 
of contacts; in the triangular lattice P(p) ~ = p 3 

for p 4. 0. We have confirmed numerically that pressure 
statistics in disordered systems are similar to the triangu- 
lar lattice; in particular, their tails have the same qualita- 
tive form and, as shown above, they obey Eq. (36 1. How- 



ever, the disordered case requires extending the theory to 
account for how pressure and tiling area are distributed 
among subpopulations with different local coordination 
numbers. This is an interesting question that we leave to 
future work. 

In light of the above discussion, we identify the single 
grain state with its pressure p and tiling area a and ap- 
proximate the entropy of the system as S = Ns, where 
the single grain entropy s is 

/>oo />oo 

s[b] = — / dp dav(p, a) [b(p, a) In b(p, a)] , (37) 
Jo Jo 



to be maximized subject to 
1 = 



dp / dav(p, a)b(p, a) 
Jo 

pOO pOO 

(p) = I dp dav(p,a)pb(p,a) 
Jo Jo 

/•OO P 00 

(a) = I dp dav(p,a) ab(p,a) . (38) 
Jo Jo 

Here v(j>, a) is the single grain density of states with pres- 
sure p and tiling area a. The result is 

b(p, a) oc exp (— ap — 7a) , (39) 

and the joint distribution P(p, a) is then 

P(p, a) = Z~ 1 v(p, a) exp (— ap — 7a) . (40) 
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FIG. 10: (solid curve) Numerically sampled average area 
(a(p)) of a tile in the triangular lattice, given that the cor- 
responding grain has pressure p. (long dashed curve) Fitted 
quadratic function 0.893((a)/(p)) 2 p 2 . (short dashed curve) 
Area of a regular hexagon with perimeter p. 



The Lagrange multipliers Z, a, and 7 are determined via 
Eq. @. 

It is convenient to factorize v(p, a) = u)(p)ip(a\p), 
where ui(p) = J dav(p, a) is the single grain density of 
states with pressure p. It is given by the righthand side 
of Eq. (35), and therefore u>(p) ~ p 3 



in the frictionless 
triangular lattice. ip(a\p) = v(p,a)/u>(p) is the density of 
single grain states with tiling area a, given that the grain 
has a pressure p. We will assume for now, and confirm 
below, that ip(a\p) is peaked at a value a*(p) « (a(p)); 
that is, given a pressure p, the most likely value of the 
tiling area (the mode of ip(a\p)) is well approximated by 
the mean area of tiles with pressure p. Under this as- 
sumption, 



da ip(a\p) exp (-7a) w exp (-j(a{p))) , (41) 



up to a prefactor that can be absorbed in Z. The 
local pressure probability distribution is P(p) — 
J daP(p, a), which becomes 



P{p) = Z V exp (-ap - 7(a(p))) 



(42) 



Thus the problem has been reduced to finding (a(p)}. 

Recall that in the reciprocal tiling, lengths are propor- 
tional to forces in the force network. In the triangular 
lattice, therefore, the pressure p on a grain is directly 
proportional to the perimeter of the corresponding tile. 
Therefore, in the simplest possible scenario, one antici- 
pates from dimensional analysis 



(a(p)) oc p 2 



(43) 



In the frictionless triangular lattice, the area a(p) of a 
tile with pressure (perimeter) p is bounded by the area 
of a regular hexagon with the same perimeter 



/ \ ^ 2 (a) 2 
a{p) -24 P =W P 



(44) 
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This is plotted in Fig. 10 which shows that the actual be- 
havior of (a(p)) is indeed quadratic, to good approxima- 
tion, and comes close to saturating the regular hexagon 
bound. From a fit to the numerically sampled (a(p)) we 
determine 



(o(p)) 



0.893-\%p 2 



(45) 



The pressure distribution P(p), determined using 
Eq. (45 1, is plotted in Fig. 11 The agreement with 
the numerically sampled distribution is excellent. Nu- 
merical distributions are computed using umbrella sam- 
pling, which allows us to determine the tail of P(p) ex- 
tremely accurately; see Appendix C for a description 
of the method. The prediction of Eq. (42 1 captures 
the cubic growth at small p (Fig. [lib), the peak near 
P ~ (p) = 6 (Fig.JTTJ;), and the Gaussian tail. The latter 
feature is best seen in Fig. Hi, which plots log P(p)/p 3 
versus p 2 . The numerically sampled distributions ap- 
proach a line with slope —1. Finite size systems fall off 
somewhat faster than a Gaussian, but deviations from 
a Gaussian tail decrease with increasing system size N. 
Note that, had the mean area of a tile (a) not been im- 
posed, we would have recovered a distribution with 7 = 0, 
i.e. an exponential tail, in clear disagreement with numer- 
ics. Thus the extra information provided by the tiling 
constraint has allowed us to capture the Gaussian tail of 
P(p). 
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FIG. 11: (color online) Numerics (solid curves) and Eq. 1 42 1 
(dashed curves) for the local pressure probability distribu- 
tion P(p) in the frictionless triangular lattice of N = 1840 
grains. The same data is shown in (a) semi-log, (b) linear- 
linear and (c) log-log plots, as well as (d) log (P/p 3 ) versus 
p 2 . To demonstrate finite size effects, (a) and (d) also contain 
data for N = 460 and N = 115. 



C. Boundaries 

All the preceding discussion has been concerned with 
the statistics of stresses in the bulk of a granular packing. 
Here we consider the influence of boundaries. 

We introduce a boundary along a triangular lattice di- 
rection and subject it to a compressive force F. In the 
orthogonal direction the packing remains periodic and 
is subject to a compressive stress such that stress ten- 
sor er is isotropic. In an experiment a packing is pre- 
pared by loading its boundary, rather than controlling 
the stress in the bulk directly. It is therefore reasonable 
to assign equal statistical weight to every microscopic 
configuration of the boundary forces consistent with the 
macroscopic load, i.e. the compressive force F. This is 
a subtle departure from the standard prescription in the 
microcanonical FNE, in which each force network con- 
sistent with a macroscopic stress carries equal statistical 
weight — all the networks of the flatly sampled FNE 
remain, but the measure in the space of force networks 
is no longer flat. A similar measure was proposed in 
Ref. [39], though there the ensemble was restricted to 
isostatic states. In the thermodynamic limit one should 
expect, and we confirm below, that the distinction be- 
tween a flat measure on all force configurations and a 
flat measure on the boundary configurations is not sig- 
nificant. We now show that, although the statistics of 



local stress in the bulk are indeed unaffected, the dis- 
tinction is important for the statistics at the boundary. 

We consider a frictionless triangular lattice with 
boundaries normal to the y-direction and periodic bound- 
ary conditions in the ^-direction. A total normal force 
Si(/b)j = F is imposed on the boundaries. For 
equally-weighted boundary force configurations {(/t>)i}, 
the distribution of boundary forces P(fb) is exponential, 
P(fb) = (/ b )- 1 exp(-/ b /(/ b )). Here (/ b ) = F/N h and 
iVb is the number of boundary contacts. Despite the 
exponential boundary force distribution, we empirically 
observe that the tail of the force distribution a short dis- 
tance from the boundary remains Gaussian, as in the 
case of a flat measure on the space of all force networks. 
Fig. [12] plots the force distribution on grains at least 5 
layers from the boundary for 8 = 2, 5, 7, and 10. The 
force distribution P(/) quickly approaches its form in the 
bulk of a periodic packing subject to a flat measure (thick 
black curve). Although the evolution of force statistics 
with depth is itself an interesting question, we leave its 
analysis to future work. 

The flat boundary force measure is interesting because 
it automatically provides an exponential force distribu- 
tion on the boundary, in good agreement with experi- 
mental boundary force measurements [4~0Ti4"2"] . Within 
the force network ensemble we find that, even with an ex- 
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FIG. 12: (a) Probability distribution P(f) of forces in the 
bulk of the frictionless triangular lattice with a boundary. 
Boundary forces are selected from an exponential distribution 
(dashed curve). Distributions are plotted for all contacts more 
than S = 2, 5, 7, and 10 layers distant from the boundary 
(thin curves). The distributions approach the form of P(f) 
in the bulk of the periodic frictionless lattice subject to a flat 
measure (thick curve), (b) Data from (a) plotted versus f 2 . 
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FIG. 13: The structure factors (a) (\p(q)\ 2 ) and (b) (\a(q)\ 2 ) 
in the triangular lattice with N = 1840 grains. The wavevec- 
tor q is parallel to a reciprocal lattice basis vector. The rescal- 
ing of the y-axis is suggested by Eq. ( 46 1 . For uncorrelated 
variables the structure factor is flat. 



ponential distribution of boundary forces, the tail of the 
force distribution in the bulk remains Gaussian. Thus, 
whether or not the flat boundary measure is realistic, it 
provides a straightforward example of a system in which 
the boundary and bulk force distributions are qualita- 
tively different. This suggests that one must be cautious 
when inferring even qualitative features of the bulk dis- 
tribution, such as the form of the tail, from experimental 
measurements on the boundary. 



D. Spatial correlations 

A single-grain picture, which was employed above to 
describe P(p) in the triangular lattice, cannot be ex- 
pected to succeed in the presence of strong interactions, 
i.e. strong spatial correlations. We therefore seek now 



to characterize the correlations in the triangular lat- 
tice. Correlations may be conveniently characterized 
by the structure factors (\p{q)\ 2 ) and (\a(q)\ 2 ), where 
ip(q) = iV -1 '}2p<f{r)e' bq ' r is the spatial Fourier transform 
of the position-dependent function (p(r). A flat struc- 
ture factor indicates the absence of spatial order. Fig. [i~3| 
shows that the pressure and area structure factors for the 
triangular lattice is nearly flat, confirming that correla- 
tions are indeed weak. 

The results of Fig. [13] may be partially motivated by 
considering the nontensile constraint. Each contact can 
sustain a normal force / which must be compressive; un- 
der our sign convention, this corresponds to a positivity 
constraint / > 0. A weaker constraint, which follows 
from positivity of the forces, is positivity of the pressure 
p on each grain. The converse is not true; pressure pos- 
itivity does not guarantee force positivity. Nevertheless, 
arguments derived from pressure positivity are useful for 
developing intuition. 

A stress state composed of the mean pressure (p) = 
Po modulated by a single oscillatory mode pp — po + 
Pq-cos (q-r) must obey \p$\ < po. Noting that the re- 
sulting constraint on fluctuations is independent of q, 
we now assume that there is a typical scale p for each 
|Pq^o| an d make a random phase approximation. Re- 
quiring ((p$-po) 2 } <Pq then gives 



Po " 



(46) 



A similar argument can be made for the area fluctuations. 
The ordinate axes in Fig. [13] have been rescaled to show 
that Eq. ((461 holds. 



Though spatial correlations are weak, they do have 
some influence on the local stress statistics. We now con- 
sider further the form of the area function (o(p)), which 



must be determined in order to predict P(p) via Eq. (42 1. 
We argued above that (a(p)) oc p 2 is to be expected both 
on dimensional grounds and because (a(p)} is bounded 
by the area of the regular hexagon with perimeter p. We 
now show that, in a system truly devoid of interactions, 
one indeed has purely quadratic scaling of (a(j>)), as in 
Eq. ( 43 ) . We also show that there are in fact small cor- 



rections to quadratic scaling, which can therefore be at- 
tributed to spatial correlations. 

Recalling that forces in a frictionless sytem are directed 
along contact normals, it is clear that the restriction to 
noncohesive forces also requires a > 0. For grain 1 with 
bond vectors {r\j} there is also a maximum possible area 
QWx(p; {fijY)- We have already noted that, in the tri- 
angular lattice, this maximum area corresponds to the 
regular hexagon with perimeter p. When interactions 
are neglected, each single grain microstate has equal a 
priori probability, independent of p and a. Moreover, 
whether the packing is ordered or disordered, scaling all 
the forces from a particular microstate of a grain by a 
scalar A > produces a new force balanced state such 
that p — > Xp and a —> X 2 a. Therefore, given the set 
of balanced states available for a particular p > 0, we 
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FIG. 14: (a) Rescaled average area p 2 {a(p)} = 
p~ 2 daa } i>(a\p) of a tile from a grain with pressure p. Di- 
viding by p 2 makes it apparent that (a(p)} interpolates be- 
tween two quadratic scalings. At small p the coefficient ap- 
proaches that predicted by a calculation neglecting spatial 
correlations. For large p the coefficient approaches the up- 
per bound given by a regular hexagon of perimeter p. (b) 
Statistics of tiles with an area a given that the tile has pres- 
sure (~ perimeter) p. For smoother curves the cumulative 
distribution C p (a) — f Q da' P(a'\p) is plotted, (thin curves) 
Cumulative distributions for p = 2,4,6, ... , 20, obtained us- 
ing umbrella sampling, (thick dashed curve) C p (a) of a single- 
grain state absent correlations with neighboring grains. Curve 
determined by numerically sampling all single grain force bal- 
anced states with a fixed pressure. Dashed vertical lines cor- 
respond to the asymptotes in (a). 



may produce all the states for p' simply by scaling each 
microstate by A = p' /p. The implication is that the con- 
ditional density of states ?p(a\p) satisfies 



iP(a\p) = \ 2 ip(\ 2 a\\p) 



(47) 



Substituting this relation in (a(p)) — da aip{a\p), one 
finds (a(p)) = cp 2 for some constant c. For packings 
under compressive stress (a) > 0, and therefore c > 0. 

Because (a(p)} must be quadratic in systems without 
interactions, deviations from quadratic scaling are evi- 
dence of interactions. To show that there are (weak) 
interactions in the triangular lattice, we calculate the co- 
efficient c directly assuming their absence. Namely, 



(o(p)> = tT 1 / d 6 /a(/i,/2,/ 3 ,/4,/5,/ 6 ) x 



(48) 



\i=l 
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where a(/i, f 2 , fo, fi, h, h) is the area of a 
tile given the grain's six forces and r\ = 

I d 6 / S (E-=i £) 8 (E Li fi-p)- Replotting the 
data of Fig. |10| by dividing out a quadratic scaling in 
p, as in Fig. 14 1, reveals that (a(p)} in fact interpolates 



FIG. 15: (color online) (a-c) Pressure distribution computed 
in the frictionless fee lattice. For a distribution decaying 
asymptotically as P(p) ~ exp (— p b ), a plot of p~ b log 10 P(p) 
will approach a flat line. Here b = 1.5 and 1.7 (see legend). 



between two quadratic scalings. For asymptotically 
small p, the area function is given by Eq. (48), while 



for asymptotically large p it obeys Eq. ( 44 1 in equality. 



Therefore Eq. ( 45 1 should be interpreted as an effective 



scaling that compromises between these two quadratic 
scalings. 

Similar behavior can also be seen in the condi- 
tional probability distribution P(a\p). For smoother 
curves, Fig. 14 3 plots the integrated function C p (a) — 
Jq da' P(a'\py If there were no spatial correlations in 
the system, plotting C p (a) against a/p 2 would collapse 
C p to a master curve independent of p. This master 
curve is the single grain density of states ip(a\p) in the 
absence of correlations, which can be obtained directly 
from Monte Carlo simulation of single grain force bal- 
anced states. Deviations from the master curve indicate 
that there are some correlations in the system, consis- 
tent with the structure factors in Fig. [13] As expected 
from consideration of Fig. |14fc ,, for small p the cumula- 
tive distribution is in good agreement with the master 
curve. For asymptotically large p the function C p (a) 
approaches a step function near the upper bound cor- 
responding to regular polygons, and C p (a) rises steeply 
over the whole range in p. This validates the approxima- 
tion that J daip(a\p)e~ 7a 



e -7<a( P )\ i.e. Eq. Un 



E. Higher dimensions 

The form of the tail of P(f) in the force network 
ensemble depends on dimensionality d, as reported in 
Refs. [HI [ID] - There the authors find P(f) decays asymp- 
totically as exp(-f bt - d ^), where b = 2(1), 1.7(1) and 
1.4(1) for d = 2, 3, and 4, respectively. The arguments 
presented above, in particular the tiling area A, are spe- 
cific to two dimensions. We now briefly consider dimen- 
sions d > 3. 

A theorem due to Minkowski states that for every set 



of d-dimensional vectors {fij}, j = 1 . 



> d, such 
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that Y]j fij = 0, there exists a unique polyhedron in d 

dimensions with the property that the unit vectors {fij} 
give the outward normals to the Zi faces of the polyhe- 
dron and the scalar magnitudes {fij} are the (d — 1)- 
dimensional areas of the corresponding faces [43] . The 
vector contact forces on each grain satisfy V . = 
absent body forces, regardless of dimension, and Z{ > d 
if grain i is not a rattler, so with each grain we can always 
associate a unique "Minkowski polyhedron" . Clearly for 
d = 2 the corresponding polyhedra are the tiles of the 
Maxwell-Cremona diagram. For d > 3 the polyhedra 
do not tesselate space; it is straightforward to construct 
counterexamples in the frictionless simple cubic lattice 
[44) . Thus the method whereby we demonstrated con- 
servation of tiling area in d = 2 cannot be generalized to 
d > 3. We stress, however, that tesselation is merely a 
means to an end; absence of a tesselation does not imply 
that the sum of polyhedra volumes cannot be conserved. 

Let us consider the consequences of a conjecture that 
in higher dimensions there exists a conserved quantity 
j{( d ) = ^2 i af\ where af^ is the d-dimensional volume 
of the Minkowski polyhedron of the contact forces on 
grain i. Then it foll ows by the entropy maximization 



arguments of Section 
f — t oo, where b(d) 



IIB2 
d/ld 



that P(p) ~ exp(-p b ( d )) as 
- 1), in reasonable agreement 



with the results of Ref. [8] . This is tested in Fig. 15 which 
plots p~ b log 10 P(p) from the frictionless fee lattice for 
both b = 1.5 and b = 1.7. Although the plot is natter over 
a broader range in p for b — 1.7, the difference between 
the two curves is small and a slow approach to a tail 
consistent with b — 1.5 cannot be ruled out. 



III. DISCUSSION AND OUTLOOK 

We have demonstrated that ensembles of hyperstatic 
force networks subject to constraints of mechanical equi- 
librium can be described within a statistical mechanics 
framework. For a given contact geometry, the ensemble 
can be explored via force rearrangements that respect 
local force balance. These rearrangements are closely re- 
lated to floppy modes, and their number is in proportion 
to the distance to isostaticity, Az. On a macroscopic 
scale, the number of force rearrangements governs the 
fluctuations of stress, and therefore pressure fluctuations 
in the FNE diverge in the isostatic limit. This divergence 
is characterized by a length scale £ w ~ Az -1 / d . 

Local stress statistics can also be explored within the 
force network ensemble, and we have extracted consid- 
erable details regarding the distribution of grain scale 
pressures, P(p). In the limit of small pressures the dis- 
tribution scales as a power law P{p) ~ p z ~ d ~ 1 , with an 
exponent that reflects the local connectivity of the net- 
work and local force balance constraints. In the limit of 
large pressures the distribution displays a Gaussian tail 
in two dimensions; this may be understood as a conse- 
quence of the invariance of the reciprocal tiling area A, 



which is quadratic in the forces. 

The force network ensemble is a minimal model; it 
is useful to the degree that it captures features of static 
granular matter in simulations and experiments, but also 
insofar as it points out necessary ingredients of more 
realistic theories. In this sense it is complementary to 
recent experimental attempts to identify relevant state 
variables in static |45j and driven [46] athermal systems. 
One striking feature of the FNE is the Gaussian decay 
(in two dimensions) of the probability density of local 
pressures p or forces /. This contrasts with early mea- 
surements of boundary forces [30-42 , though we saw in 
Section II C| that statistics at the boundary need not cor- 
respond to bulk statistics. Early theoretical efforts such 
as the q-model also predicted exponential tails but only 
incorporated scalar force balance [47] ; we have seen that 
Gaussian tails are intimately connected to the reciprocal 
tiling, which requires vector force balance. We consider 
the form of the tail of P(f) in real granular systems an 
interesting and open question; Ref. [IU] summarizes re- 
cent theoretical, numerical, and experimental work on 
the matter. 

The FNE is easily extended to frictional packings, and 
this is an obvious avenue of future research. The ensem- 
ble may also be used to study departures from equal a 
priori sampling of states, which is permitted - even likely 
- in non-equilibrium ensembles. 
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Appendix A: Degree of freedom counting 

Conventionally, the microcanonical FNE on a given 
frictionless sphere packing has been defined as the set 
of force networks satisfying the matrix equation 



Af = b, 



(Al) 



and a set of inequalities [3J 03 13 H5J. The inequalities 
restrict f to noncohesive states in which each fij > 0. 
A is a zN/2 x (dN + d(d + l)/2) matrix. Its first 
dN rows encode force balance on the given contact net- 
work, i.e. y)j Aij\f]j = Fix, the x-component of the 
net force on grain 1, £\ ^2j[f]j = Fly: an d so on - The 
final d(d + l)/2 rows of A are chosen so as to return 
each of the unique stress tensor components of f , namely 
o a p = (l/2V)X^[^aL/y,3- Correspondingly, the first 
dN elements of b are zero: each grain experiences no 
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net force (here we assume for simplicity a periodic pack- 
ing absent body forces). The final entries of b contain 
the desired values of the stress tensor components. The 
matrix A is underdetermined, and therefore solutions to 
Eq. ( Al ) are of the form 



fo 



f A 



(A2) 



fo is a particular solution to Eq. (Al) and the {f„} are 
null ve ctors of A. Aa is the nullity of A. The similarity 
of Eq. (A2) to Eq. (10) is obvious, and in fact what we 



have done in Section |I| is to construct the null vectors of 
A "by hand" , without writing down A itself, simply by 
considering the constraints of local force balance - the 
force rearrangements {£f° xt } preserve local force balance 
and the stress tensor components, and therefore are null 
vectors of A. Here we show that = N w , i.e. that 
the number of null vectors of A does indeed equal the 
number of force rearrangements we identified. 

In a periodic static packing, zero net force on the entire 
system follows directly from periodicity. Hence J^. F t = 

T,i Y,ifij =0, or 



N-l 



(A3) 
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-/c+l 
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FIG. 16: (a) Indexing of the Zi contacts around grain i sub- 
ject to frictional contact forces (black arrows) . (b) Reciprocal 
tile of the grain in (a). Coordinates of the tiling vertices are 
labeled {h c }. Vertex h c is connects edges c and c+l (modulo 
Zi). (c) The building block of a frictional force rearrangement 
involves four grains a, b, c and d, and two voids 1 and 2. (d) 
Coordinates of the tiling vertices hi and /12 corresponding to 
voids 1 and 2 of (c). 



one of the w's can be expressed as a linear combination 
of the other N-l. 



Therefore it suffices to impose force balance on N — 1 of 
the grains; force balance on the final follows "for free". 
This fact can be used to count the packing's (force) de- 
grees of freedom. A periodic triangulation with N ver- 
tices has 3A edges, and hence supports 3N compressive 
forces. There are 2(A — 1) force balance constraints in a 
two-dimensional system, hence 3N — 2(N — 1) = N + 2 
degrees of freedom in the forces. Further constraining 
the three components of the stress tensor (a microcanon- 
ical FNE) leaves N\ = N — 1 degrees of freedom. If 
the contact network is not a triangulation there is one 
less degree of freedom for each deleted edge, namely 
N A = N - N d - 1 = (l/2)AzN- 1. (Recall z = 6 in 
a triangulation.) Therefore Aa = N w . Though we have 
not proven linear independence of the disordered wheel 
moves, we have checked for, and always found, linear in- 
dependence in numerically generated triangulations. 

A w = (1/2) Az N- 1, rather than (1/2) Az A, because 
there is a null direction associated with the N wheel 
moves in a triangulation. We now construct that null 
direction. In a triangulation each vertex i has an associ- 
ated local pressure difference 6pi, expressed with respect 
to the pressure on grain i in reference state fo- It is 
straightforward to see that the {Spi} can be formulated 
as a linear superposition of the disordered wheel move 
weights, Spi — LijWj. (In fact the square matrix 
is a discrete representation of the Laplacian operator.) 
There are N <5p's, but one is redundant because in the 
microcanonical FNE Spi = for every force network 



Appendix B: Rearrangements in frictional packings 

In frictional disk packings it is again possible to con- 
struct force rearrangements from localized objects in the 
Delaunay triangulation. Although these objects are no 
longer floppy modes in the reciprocal tiling, the tiling is 
again helpful in identifying them. Once the local rear- 
rangements have been identified, they can be combined 
to produce spatially extended force rearrangements in the 
packing in a manner completely analagous to Section [I E| 
As in the frictionless case, we will identify a set of 
motions of the tiling vertices that constitute local rear- 
rangements of the forces. Because the motions are in 
the tiling, the rearrangements are guaranteed to respect 
force balance. In the frictionless case, the additional con- 
cern was ensuring that, after rearrangement, all forces 
remained normal to their contacts. In the frictional case 
this restriction is lifted, but we must now be sure that 
the rearrangement respects local torque balance. 

Consider a disk as in Fig. [l6fc , with contacts located 
at {R c } and forces {/ c } acting on those contacts. With- 
out loss of generality we place the origin at the disk's 
center. Torque balance requires ^ c R c x f c = 0. Follow- 
ing Ref. [28], it is convenient to rewrite this expression 
in the following way. Label the z contacts of the grain 
c = 1 ... z in a righthand fashion (Fig. 
bel the z vertices of the grain's reciproca 
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in the ensemble (see Section IF). Therefore the {wj } sat 



isfy the sum rule lj 



0, where lj — '^2 i Lij, and 



Also la- 
tile c = 1 . . . z 

such that contact c connects vertices c and c+l (modulo 
z). We first construct the set of vectors {g c } depicted 
in Fig. 16 1 (dashed arrows), where g c := R c +i — R c . 
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Secondly, note that the contact force f c is the difference 
of tw o ve rtex coordinates h in the tiling rotated by ir/2 
(Fig. 16 d), i.e. f c = —(h^ +1 — hjr). Using these defini- 
tions, the torque balance constraint on the grain is 



xhir = 0. 



(Bl) 




(a) 




Eq. (Bl) is useful in constructing a localized force re- 
arrangement. Consider the configuration of four grains 
labeled a to d in Fig. 16 There are six relevant g vec- 
tors. Adjusting our previous notation, we label, e.g., gia 
the vector around void 1 contained in grain a. Similarly, 
there is a tiling vertex associated with each of the voids; 
we label them hi and h 2 . We wish to identify a change in 
each, Shi and Sh 2 , that respects torque balance. Such a 
change will enter the torque balance constraint, Eq. (Bl I, 
on each grain in the figure. There are then four con- 
straints, 

= Qla x Shi 

= Qu x Shi + g 2b x 5h 2 

= Q 2 c x Shi 

= Qid x Shi + Q2d x Sh 2 ■ (B2) 

Hence there are four unknowns, the components of Shi 
and 5h 2l and four constraints. As in the frictionless case, 
however, there is a degeneracy in the constraints. Con- 
sider the sum of the four constraints of Eq. (B2). One 
finds 

= (Qla + Qlb + Qld) X Shi + (92b + Q2c + Q2d) x 8h 2 ■ 

. (B3) 

This equality holds independently of Shi and Sh 2 be- 
cause the terms in parentheses are sums around closed 
loops. Therefore Eq. (B2) represents only three indepen- 



FIG. 17: The counterparts to wheel moves in (a) the frictional 
triangular lattice and (b) the frictional square lattice. Arrows 
indicate the vector changes in force at each contact. 



dent constraints, and the resulting one parameter family 
of motions is a localized force rearrangement in a fric- 
tional packing. 

Every edge in a Delaunay triangulation is shared by 
a pair of triangles, hence we can identify one such re- 
arrangement for each of the Az N edges in the triangu- 
lation. For general disk packings, these frictional dis- 
ordered wheel moves can be constructed in the pack- 
ing's Delaunay triangulation and then combined to form 
extended force rearrangements that change neither the 
normal nor the tangential force components on deleted 
edges; this is done in a manner entirely analogous to the 
treatment in Section |U 

In Fig. [17] we depict the equivalent of wheel moves in 
two ordered frictional contacts networks, the triangular 
and square lattices, which were employed in Ref. [9]. The 
move in the triangular lattice corresponds directly to the 
situation in Fig. |16( while the move in the square lattice 
is a superposition of multiple elementary building blocks. 




FIG. 18: (a,b) Probability distributions of the maximum force 
/max in a force network f and (c,d) contact force distributions 
P(f) for a frictionless triangular lattice (TV = 1840). Three 
different situations are considered: (i) the force network en- 
semble without umbrella sampling (fne), (ii) the ensemble n in 
which W(f m ax) is chosen such that P^(/ ma x) is approximately 
flat (it), and (in) the ensemble n in which ensemble averages 
are reweighted to the force network ensemble using Eq. |C4| 
(fne-7r). In all cases, the results for umbrella sampling (fne-7r) 
are identical to those computed without umbrella sampling 
(fne). All forces are normalized such that (/) = 1. Con- 
tact forces larger than 5 (/} are hardly sampled in the force 
network ensemble unless umbrella sampling is applied. 



Appendix C: Umbrella sampling 

In simulations we employ umbrella sampling, which 
permits extremely precise determination of the probabil- 
ity density of large stresses [8] . Monte Carlo simulations 
sample force networks f with a probability proportional 
to their statistical weight, i.e. in the force network ensem- 
ble each force network that satisfies local force balance 
and the boundary conditions is equally likely. As the vast 
majority of force networks only contain contact forces of 
order of magnitude (/), large contact forces or local pres- 
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FIG. 19: (a) Typical force network containing a large contact 
force. Edge thicknesses are proportional to force magnitudes, 
(b) Reciprocal tiling of the force network in (a), (c) Typical 
force network containing a large local pressure (d) and its 
tiling. Colors map redundantly to force magnitudes or tile 



sures are hardly sampled and it is not possible to obtain 
P(f) or P(p) accurately for large / or p respectively. 

To improve the sampling for large contact forces or lo- 
cal pressure, we employ the umbrella sampling method 
[48| [49] . The central idea is to create a bias in the force 
networks obtained in the Monte Carlo simulations, and to 
correct exactly for this bias afterwards. To illustrate this, 
let us denote the a priori probability of a force network 
f by G(f). In the force network ensemble G(f) equals ei- 
ther or 1. Monte Carlo simulations of the force network 
ensemble generate configurations with a probability pro- 
portional to G(f) so therefore the average of a property 
A can be computed from 



{A) = 



K 



(CI) 



in which fi, f2, - • ■ , (k are the force networks generated 
by the Monte Carlo scheme. To generate more force net- 
works with large forces, consider the ensemble tt in which 
the a priori probability of a force network f equals 



G 7r (f)=G(f)exp[W(f)] 



(C2) 



in which W(f) is an arbitrary function that only depends 
on the force network f . Monte Carlo trial moves from 
state fo to state f„ in this ensemble are thus accepted 
with a probability [50 



acc(fo — ¥ f„) = min I 1 



G(to) 



exp[W(f n ) - W(i )} 



Ensemble averages calculated in the ensemble tt can eas- 
ily be converted to ensemble averages in the original force 
network ensemble, as 



(A) = 



JdiG(i)A(i) 
JdtG(f) 

J dfG(f)exp[W(f)]A(f) exp[-W(f)] 



/ dfG(f) exp[W(t)] exp[ 
(A(f)eM-W(i))^ 
(exp[-W(f)]) w ' 



-W(f)] 



(C4) 



in which we used the shorthand (■••)_ for averages com- 
puted in the ensemble tt. A smart choice of W(f) will 
sample many networks with large contact forces, so P(f) 
can be computed accurately for large /. A convenient 
choice is to introduce the order parameter / max (f) as 
the largest contact force of a force network f and to 
express W as a function of / max only. The function 
W(/max(f)) can be chosen such that the probability dis- 
tribution P 7r (/ m ax) (computed in the modified ensemble) 
is approximately flat. As from Eq. C4 follows that 



P(f max) = constant x P^(f max ) exp[-W (f max )] , (C5) 



it is convenient to iteratively determine lf(/ max ) such 
that 



W(f n 



lnP(/ max ). 



(C6) 



(C3) 



To illustrate the use of umbrella simulations, in Fig. [18] 
we have plotted the probability distributions of P(f) and 
P(fmax) computed with and without umbrella sampling 
for a frictionless triangular lattice of N = 1840 particles. 
From this figure it becomes clear that umbrella sampling 
increases the accuracy of P(f) for large contact forces by 
many orders of magnitude. 

To compute the propability distribution of local pres- 
sures P(p) accurately for large p, we employ the same 
scheme. However, it turns out configurations contain- 
ing a single large local pressure are qualitatively differ- 
ent from configurations containing a single large con- 
tact force. In the former, the large pressure is spread 
out over all contacts of the particle, while in the latter 
only two grains experience a large contact force. This 
is shown in Fig. [l9j As a consequence, to compute 
P(p) accurately for large p a different order parameter is 
needed. In our simulations to compute P(p), we choose 
W = W(p max (f)) in which p max (f) is the maximum local 
pressure of a force network f . 
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We exclude using networks on the boundary of the poly- 
tope for f*o, which can cause difficulties with sampling. 
See Ref. [4] for a discussion of this point. 
Other sampling methods are possible. E.g. when employ- 
ing umbrella sampling (see Appendix C) moves are ac- 
cepted/rejected according to Eq. (C3l. 



Of course, the limiting case itself is uninteresting in the 
force network ensemble: for truly isostatic packings the 
microcanonical FNE contains a single force network, and 
the canonical FNE is composed of rescalings thereof. 



